Spatial epidemiology of hemorrhagic disease in Illinois wild white-tailed deer

Epizootic hemorrhagic disease (EHD) and bluetongue (BT) are vector-borne viral diseases that affect wild and domestic ruminants. Clinical signs of EHD and BT are similar; thus, the syndrome is referred to as hemorrhagic disease (HD). Syndromic surveillance and virus detection in North America reveal a northern expansion of HD. High mortalities at northern latitudes suggest recent incursions of HD viruses into northern geographic areas. We evaluated the occurrence of HD in wild Illinois white-tailed deer from 1982 to 2019. Our retrospective space–time analysis identified high-rate clusters of HD cases from 2006 to 2019. The pattern of northward expansion indicates changes in virus-host-vector interactions. Serological evidence from harvested deer revealed prior infection with BTV. However, BTV was not detected from virus isolation in dead deer sampled during outbreaks. Our findings suggest the value of capturing the precise geographic location of outbreaks, the importance of virus isolation to confirm the cause of an outbreak, and the importance of expanding HD surveillance to hunter-harvested wild white-tailed deer. Similarly, it assists in predicting future outbreaks, allowing for targeted disease and vector surveillance, helping wildlife agencies communicate with the public the cause of mortality events and viral hemorrhagic disease outcomes at local and regional scales.

Other wild ruminant species such as mule deer (Odocoileus hemionus), pronghorn (Antilocapra americana), North American elk (Cervus canadensis), and bighorn sheep (Ovis canadensis) can contract EHDV; however, WTD are particularly susceptible to severe disease associated with EHDV 26,27 . During a severe EHD outbreak in captive WTD in Missouri in 2012, the average within-herd morbidity rate was 46%, while within-herd mortality rates ranged from 3 to 84%, with an average of 42% 28 . Herd mortality in wild WTD in Missouri was estimated at up to 16% 29 ; however, in locations where deer populations may not be regularly exposed to EHDV, mortality can be much higher 30,31 . Because of the potential for large die-offs during severe outbreak years, EHD remains an important viral disease of WTD in Illinois and North America 32 .
Population-level studies of HD on wild WTD would assist in the development of strategies to help manage a public resource that impacts local economies, especially because HD outbreaks often precede or coincide with hunting seasons. Despite recognizing the disease since the late 1800s in the USA and increasing reported cases in recent decades, HD's spatial and temporal distribution in wild WTD is understudied. Spatial epidemiology methods have been used to analyze infectious diseases in wildlife, including chronic wasting disease in Illinois WTD 33 , HD in WTD in Missouri 34 , and Baylisascaris procyonis (raccoon roundworm) in southern Ontario, Canada 35 . Analyzing large datasets using retrospective space-time permutation models helps focus surveillance efforts on high-risk areas to detect disease outbreaks 36 .
Past HD studies in the Midwest focus on region-wide trends in HD expansion, the detection and occurrence of new HD serotypes in the USA, and specific outbreak years in domestic and wild ruminants 5,6,28,37 . However, studies on HD in wild ruminants in the Midwest are still lacking 8,11 . Recent studies suggest HD outbreaks are occurring in more northern latitudes 5,38 . However, most HD data in wild WTD often report only the annual, county-level presence or absence of disease 8 rather than numerical counts of dead deer. Thus, characterizing outbreaks in wild herds across time and space is difficult. Space and time clusters of local HD outbreaks in wild WTD populations could reveal expanding virus and vector distribution in Illinois that could be caused by changes in land use (e.g., habitat loss or agricultural development) and climate. Spatial and temporal modeling of HD outbreaks in wild WTD would highlight research gaps in the epidemiology and surveillance of HD in Illinois and identify future research directions.
BTV infection in cattle and other domestic livestock is well studied in the USA, and there is evidence of BTV seroprevalence in Illinois cattle 39 . However, evidence suggests that BTV infection in Illinois WTD may be low 40 . In contrast, the three EHDV serotypes (EHDV-1, -2, -6) known to occur in the USA have been documented in Illinois WTD 37 . The once exotic EHDV-6 was first isolated from wild and captive WTD in Indiana and Illinois in 2006, followed by Missouri, Texas, and Kansas between 2006 and 2008 37 . During fall of 2013, the first EHD outbreak in domestic ruminants was reported in Illinois. At least 18 cattle herds in northwestern Illinois were affected 41 . This outbreak coincided with the EHD outbreak affecting wild WTD in that region 41 and followed the major outbreaks in domestic ruminants in the USA in 2012 28 . Cattle and deer have been considered important vertebrate hosts in the EHDV/BTV transmission cycle 9,42 . Furthermore, recent studies in the USA have demonstrated animal host density 43 and biting midges preference patterns for different Cervidae and Bovidae hosts 44 as important factors for maintaining HD viruses in captive ruminant herds. The occurrence of HD in captive herds may have consequences for subsequent infection in wild WTD populations 9 . In Illinois, cattle and deer often inhabit the same large-scale geographic areas and are often close to habitats where HD competent vectors occur. However, because of the severe impact of HD and the unknown occurrence of BT in wild WTD in Illinois, targeted surveillance efforts are needed to identify risk factors contributing to HD outbreaks in wild populations.
Currently, there are no studies on the historical occurrence of EHD or BT in wild Illinois WTD. Spatial epidemiology can identify areas with increased disease rates. Therefore, a preliminary investigation of increased HD rates to identify areas and periods with higher-than-expected HD reports will serve as a baseline for future studies. Additionally, such baselines can guide future studies to evaluate additional gaps in knowledge related to local risk factors, predominant serotypes, competent vectors, and potential animal-vector-environmental interactions favoring HD outbreaks in wild WTD populations. Therefore, our objectives are (1) determine the distribution of HD reports in wild WTD by county in Illinois and characterize the intensity and temporal intervals of HD outbreaks, (2) identify records that indicate the presence of BTV in wild Illinois WTD, and (3) evaluate geographic expansion trends of HD in Illinois using a space-time scan statistic model. To achieve the analytical components of the study objectives, we constructed a choropleth map illustrating county-level numbers of HD cases and an isopleth map to account for the areal bias associated with county boundaries. Furthermore, we used a hot-spot analysis to identify counties with significantly higher and lower numbers of HD cases. Finally, we used a space-time permutation scan that offers a time component for the identified HD clusters in addition to the spatial statistical steps.

Results
Descriptive analysis. From 1982 to 2019, 99/102 Illinois counties reported HD at least once, and three counties (Boone, DeKalb, and Kankakee) had zero reports of HD. One county (Grundy) reported presence of HD but not the number of deer with HD cases (Fig. 1). No HD cases were reported in 14 of the 38 years studied (1982-1987, 1990-1991, 1993-1995, 1997, and 1999-2000). We documented 8044 total HD deer (including reported dead deer from outbreaks and confirmed HD deer from diagnostic results). There was an average of 4. 6  Disease mapping. The total number of HD cases/county for the full study period ranged from zero to 523 cases. Counties with a high number of HD cases were in west-central, southeast, and northeast Illinois (Fig. 3a). Empirical Bayesian kriging results are illustrated in an isopleth map (Fig. 3b). Several regions with an increased number of HD cases were detected, including the west-central, northeast, and south regions of Illinois.
The hot spot analysis identified five counties with a low number of HD cases (cold spots) in north-central Illinois and nine counties with a high number of HD cases (hot spots) in east-central Illinois (Fig. 3c).
Space-time scan statistic. The retrospective space-time scan statistic detected four high-rate clusters ( Fig. 3d; Table 1

Discussion
This spatial epidemiology study allows identifying high-risk areas supporting future multidisciplinary research that will help determine potential animal-vector-environmental interactions during outbreaks of HD on wild WTD populations. Furthermore, it will guide field studies to evaluate the serotypes, competent Culicoides spp., and host reservoirs present at local scales.     (Fig. 3a). These findings of northern expansion are consistent with previous reports 5 and with the high morbidity and mortality usually reported when naïve populations are exposed to HD viruses 9 .
During the timeframe of each cluster, Illinois experienced drought in the late summer and fall months, especially in 2012 when the state experienced extreme drought conditions along with much of the US 5,28,46,47 .
Although not evaluated in this study, the correlation of HD outbreaks with drought and other climatic variables (e.g., temperature, humidity, wind) has been a focus of recent studies 5,21,23,48,49 . Drought has been suggested to contribute to HD outbreaks 21,28 . Warm and shallow water with freshly exposed mud is ideal for midge reproduction, especially C. sonorensis 50 . The congregation of animals due to limited water sources may contribute to increased exposure to Culicoides 50,51 , mud disturbance, and fecal matter deposition, increasing habitat suitability for propagation and reproduction of Culicoides vectors 9 . The role that drought and other climatic factors may play in vector dynamics, virogenesis, and HD outbreak severity could be important considerations in anticipated changes in climate and hydrology 52,53 . Although the clusters identified in this study were associated with years where drought was reported in Illinois, not all HD outbreaks in the USA are related to drought conditions, as demonstrated in the eastern USA in 2017 54 . Therefore, along with climate impacts on vector activity, other biotic and abiotic factors may impact HD transmission and outbreak severity. Landscape characteristics 55 , soil chemistry 17,56 , land-use practices, and environmental conditions associated with livestock 18,57-60 (e.g., overflow from stock tanks, livestock wastewater lagoons, soil disturbance along riparian areas or stock ponds) may also impact vector abundance and vector-host interactions. A better understanding of how these and other factors affect the disease system could inform risk communication and management decisions to guide hunter harvest in outbreak areas and sustain healthy populations of wild deer in severe outbreaks.
Many HD positive counties are located along the Illinois River and are included in three of the four clusters (C2, C3, and C4) identified by the space-time analysis (Fig. 3d). Pike County was identified as a cluster (C4) in 2006 and is part of the nine counties identified as hot-spots for HD cases in Illinois. This county reported HD in nine of the 38 years of data. Pike County is a mix of forested and agricultural land between the Mississippi and Illinois Rivers (Supplementary Video 1). This county consistently reports some of the highest annual deer harvest numbers in the state 61 , which may be due to higher deer densities, increased hunting participation, and the perceived "trophy quality" of Pike County WTD. These same factors may also contribute to higher reporting of HD because the public may be more invested in this resource, thus contributing to increased surveillance and sampling efforts. Higher deer populations may also contribute to the congregation of reservoir ungulate hosts, facilitating exposure to infected Culicoides 62 . Because of higher deer populations 61 and possibly because of the availability of vector habitat (e.g., wetlands, stock ponds, and exposed riparian areas), Pike County and counties with similar physiographic characteristics along the Illinois River may be more prone to HD outbreaks than other parts of the state. The study of animal movement around water sources and higher coverage of water bodies have not been evaluated in association with HD in Illinois. While this study did not evaluate river corridors and their relation to HD outbreaks, we recognize that defining the role of animal movement around water sources and the link of the amount and type of water bodies by county on disease spread may help determine inter-related variables affecting HD dynamics. Incorporating these variables in future HD epidemiological models may contribute to understanding the changing patterns and distribution of HD in Illinois.
Previous studies have identified multiple factors affecting host-vector interactions and vector distribution and abundance. High animal densities in captive populations-independently if animals are kept in a fenced enclosure, open pasture, or free-roaming in natural preserves-have been identified as significant contributors to HD presence compared to wild populations in the same geographic area 43 . Additionally, the preference patterns of different Culicoides for different ungulate species have also been suggested as a contributing factor in HD dynamics between ruminant hosts 44 . Facilities with exotic ungulates, although not developing the disease, could contribute to local EHDV/BTV transmission 43 . The high density of animals can impact landscape characteristics and contribute to higher vector density 63,64 . Other factors such as climate and weather can affect vector ecology. Although unclear which physiographic features could be responsible for the increase in reported HD cases during the outbreaks in Illinois, identifying the hot-spot/cold-spots in specific areas and during a period can help www.nature.com/scientificreports/ guide future surveillance and management efforts at the county level. Further investigation of climatic features, Culicoides spp., artificial changes to the environment (e.g., urban development and agricultural practices), animal movement between farms, and other risk factors in hot-spot areas can help identify potential predictors for future outbreaks in specific areas in Illinois. Recognizing the specific risk factors by area will help develop more specific management plans to control disease dissemination in domestic ruminants, thus protecting other domestic and wild ruminant populations. While we did not evaluate geographical factors affecting HD reporting, deer that die from HD may be more likely to be detected in areas closer to urban settings than those in rural settings where they may go undetected. Host distribution differs between urban and rural areas, and therefore, die-offs are more likely to be reported in areas with deer habitat and greater host densities. Outdoor activities such as farming, hiking, or hunting may increase HD reports. Hemorrhagic disease data in Illinois are obtained through diagnostic submissions, the Southeastern Wildlife Cooperative Disease Study (SCWDS) annual HD questionnaire, and available harvested deer samples tested for HD outside of the outbreak season. Historical data pertaining specifically to HD in Illinois wild WTD was difficult to obtain. The earliest available reports of HD in Illinois date to 1982. Since 1982, SCWDS has conducted a nationwide HD survey of state wildlife management agencies. The goal of the survey is to monitor county-level reports of suspected and confirmed HD in wild WTD and other wild ruminants. However, the survey does not consistently capture the number of deer affected in each county or state, but it may capture county-level presence. The Illinois Department of Natural Resources (IDNR) provided annual HD reports beginning in 1998, which included county-level counts of reported dead deer. Beginning in 2013, available IDNR data included more detailed information regarding case reports (e.g., location, date, sex, and age), although this is not always included in the records. We recognize the challenges and limitations to conducting HD surveillance during outbreaks, which include: (1) the time of year (late summer to early fall) when decomposition may accelerate due to high temperatures, (2) location of carcasses (commonly found in water), and (3) limited availability of field staff for sampling a high number of dead deer during large-scale outbreaks. All these factors make it extremely difficult to collect the necessary tissue samples for HD diagnostics in a timely manner (within 24-48 h of a death to avoid decomposition), impacting sampling efforts over the years.
There is insufficient information to know how much the sampling effort varied over the years. While recognizing the challenges and difficulty in gathering this information for every HD dead animal during an outbreak, variations in the data make complex data analyses challenging. To advance the quantification and analysis of risk for an outbreak of HD in wild WTD and to permit the development of better predictive models for HD outbreaks, officials should standardize the data collected from HD reports and field visits to include: date of reported case, county, TRS (township, range, and section number), geographic coordinates, number of dead deer reported, age (fawn, yearling, or adult), sex, and sample ID (if tissues are submitted to a diagnostic laboratory).
Using available Illinois historical HD data, we identified diagnostic results for wild WTD in Piatt County, Illinois, with antibodies for BTV. Eleven deer samples from 2005 to 2006 tested positive for antibodies for BTV. Of those, five tested positive for antibodies for both EHDV and BTV. Cross-reaction is known to occur when using the AGID test to identify BTV and EHDV 65 . However, because six of those deer samples tested positive for BTV antibodies and negative for EHDV antibodies, cross-reaction may not have occurred in these cases. The deer that tested positive for antibodies against BTV were all harvested from the same area in Piatt County. These Piatt County deer and two additional deer found positive for BTV in Peoria and St. Clair counties from a previous study 40 demonstrate that BTV exposure has likely occurred in wild WTD in different regions of central Illinois. However, available diagnostic results of wild WTD statewide suggest that EHDV dominates HD mortality in Illinois deer. None of the deer tissue samples submitted during the outbreak months tested positive for BTV. This agrees with findings in domestic ruminants during the 2013 outbreak, where four northern Illinois counties (Kane, Jo Daviess, Henry, and Stephenson) reported EHDV infection in cattle, but no BTV cases 41 .Thus, while we know that BTV is present in Illinois as evidenced by seropositivity in Illinois cattle 39 and that deer become infected with BTV, we do not know its distribution or the extent to which wild WTD are infected and affected by this virus. Domestic cattle are considered an important host for BTV and EHDV 66 . However, the role of wild WTD as an important host has been associated only with EHDV 9 . However, cattle may not be the only farmed ruminant that plays a crucial role as a host for BTV/EHDV. For example, a recent study identifies 44% of mortality cases in farmed WTD in Florida between 2012 and 2020 associated with BTV, EHDV, or co-infection of BTV/EHDV 67 . This study demonstrates the impact of these viral diseases in captive cervids and highlights the importance of continued surveillance and using farmed ruminants as sentinels for HD in the surrounding wild ruminant populations 67 .
Nevertheless, our results corroborate that, despite detecting some cases with positive EHDV and BTV results, EHDV was the orbivirus driving the HD cases in wild Illinois WTD (Supplementary Table S1). Previous studies have identified that increased temperatures can affect vector survival and competence for different serotypes. For instance, temperatures between 27 and 30 °C increase vector competence for EHDV (serotype 1); however, these temperatures interfere with vector competence for BTV (serotypes 10 and 16) 23 . Although we did not evaluate the effect of temperature on the HD outbreaks in Illinois, we recognize the importance of expanding HD surveillance to identify the association of environmental variables with distribution and abundance of competent Culicoides. Expanded surveillance will help detect the incursion of new virus strains to Illinois and to North America, while establishing the most common serotypes in Illinois. Only through surveillance can a complete picture of the variables associated with the EHDV/BTV triad (vector-host-environment) be developed and used to assess management programs designed to protect wild and domestic ruminants.
We demonstrated that HD cases in Illinois cluster in space and time, suggesting that we are dealing with the localized transmission of the virus in a period. Additionally, we found an expansion of HD to areas of the state that had previously not reported cases. However, we did not evaluate ecological factors (e.g., changes in deer behavior and population dynamics and Culicoides population dynamics) and anthropogenic factors (e.g.,

Scientific Reports
| (2022) 12:6888 | https://doi.org/10.1038/s41598-022-10694-y www.nature.com/scientificreports/ landscape changes, introduction or translocation of wildlife and livestock) that may impact severe HD outbreaks. Notably, cases in earlier years were reported from central and southern Illinois, and in later years, they were also reported in northern Illinois (Fig. 1). Regional patterns in HD outbreaks have historically fallen into one of two categories: enzootic (occurring every 2 to 3 years) or epizootic (occurring every 8 to 10 years) 25,68 . Figure 1 and a time series plot of the number of counties reporting HD occurrence between 1988 and 2019 ( Supplementary  Fig. S2a) revealed three peaks of HD in 1998, 2007, and 2012. The peaks also demonstrate an increase in counties with a HD report. Visually inspecting Fig. 2, we observe a cyclical pattern with a high number of HD cases reported in 2007 and 2012-2013. This pattern is also observed in a time series plot (Supplementary Fig. S2b). This observation was confirmed by the scan statistical analysis where three large space-time clusters of the number of cases were detected in the south, north, and central regions of Illinois during these periods (Fig. 3d). Our study identified BTV and EHDV exposure in deer through general disease surveillance outside typical HD outbreak periods and highlights a potential gap. By limiting surveillance to outbreaks, we do not detect areas where BTV and/or EHDV infected deer survived. While HD outbreaks may be seasonal, surveillance efforts, particularly serological surveys, should be conducted year-round. Although we did not evaluate environmental factors and the biological vectors influencing HD transmission in Illinois, we recognize the importance of climate change on emerging and re-emerging viral diseases at a local and global scale 9,69 . For instance, the identification of EHDV-6 in Indiana and Illinois in 2006, Missouri in 2007, Kansas and Texas in 2008, and Michigan in 2009 were suggested as a result of overwintering maintenance and, therefore, expansion of EHDV in the USA 70 . However, there is no evidence of the virus overwintering in northern areas; therefore, this remains to be evaluated. During the 2006 EHDV outbreak in Israel, winds were considered a major contributing factor to HD spread 71 . Additionally, the maintenance of BTV in a farm in temperate regions in northern California was associated with prolonged survival of C. sonorensis, thus contributing to interseason maintenance 72 . Therefore, there are local and global scale indicators of environmental factors affecting HD transmission that should not be overlooked. Moreover, phylogenetic analysis of virus strains may help identify the movement of virus populations 70 and the genetic and spatial distance of virus strains identified during outbreaks. Such findings may point to other variables contributing to virus transmission, dissemination, vector competence and maintenance. HD surveillance on farmed ruminants may serve as a sentinel strategy to identify risk factors for disease transmission in captive and wild populations 67 . Livestock and captive cervid facilities may be valuable sentinels for HD surveillance because it is possible to identify the recently infected (viremic) animals necessary for the infection of biological vectors 43,67 . Evaluation of native Culicoides species near captive facilities may help detect shifts in HD occurrence and the potential for maintenance of HD viruses outside the known vector home range 9 .
In regions like Illinois, HD outbreaks may result in high mortality and impact local economies and tourism, especially the outdoor recreation industry. However, HD outbreaks also affect the agricultural industry and captive cervid herds. Recognizing the importance of HD, IDNR collects detailed information when investigating reports of dead deer. However, it is not always possible to collect geographic coordinates for each dead deer. Nevertheless, geographic coordinates are critical for fine-scale analyses of outbreaks and future studies on vector, landscape, and natural and agricultural conditions that may play a role in HD outbreaks. When the collection of tissue samples is not possible (e.g., poor post-mortem condition of the carcass) or field visits are not feasible (e.g., too many reports or limited staff), location data should be collected. TRS location data provides the next best option and offers an opportunity to evaluate problem areas where repeated outbreaks occur. Additionally, TRS-level data can assist with predictive models of disease expansion while accounting for variables like wetlands, forests, soils, and proximity to livestock facilities. Predictive models can help wildlife agencies prepare resources for outbreak response, inform the public of expected outcomes, and help inform management decisions (e.g., adjusting harvest quotas) to sustain healthy populations of deer in areas of severe HD outbreaks.

Methods
Study area. Illinois is a state in the Midwest region of the USA and is made up of 102 counties (Supplementary Fig. S1). It is roughly 150,000 km 2 , ranges in latitudes from 42.5° N to 37° N, and is approximately 630 km long from north to south. Illinois' climate is temperate and includes four seasons. The average annual precipitation ranges from ~ 1200 mm in the south and ~ 812 mm in the north. Illinois is predominantly agricultural land mixed with moderately to highly dense urban areas and some deciduous and deciduous-mixed forest land, prairie, and wetlands. HD data collection in Illinois. Currently, in Illinois, HD surveillance is achieved through IDNR investigation of public reports of dead deer. Deer reported as having died from HD met at least one of the three mortality criteria suggesting acute HD: (1) high and acute mortality events during late summer and early fall, often near a water source (e.g., river, creek, pond), (2) HD diagnosis from necropsy findings, or (3) dead deer with laboratory confirmation of hemorrhagic disease (EHD or BT) 32 . IDNR staff record information from all HD reports including the number of dead deer and the location of the outbreak (at least the county level). Field investigations are conducted whenever possible but are limited by staff availability and the condition of the carcasses reported. Priority is given to counties in which HD has yet to be confirmed in a given year. Once HD is confirmed in a county, it is generally not revisited for laboratory confirmation of other dead deer during the same season. Circumstances that affect the collection of useful diagnostic samples include temperature, location of the carcass, and age of the carcass. The level of tissue decomposition increases with time, environmental temperature, humidity (in air and from carcasses found in water) 75 . The number of days post-mortem before the carcass can be examined by biologists affects the quality of a sample available for diagnostic testing. Thus, IDNR staff collect spleen and/or lung samples from fresh deer carcasses (usually within ~ 48 h of death) for laboratory confirmation (IDNR personal communications). Reported HD cases are logged by IDNR and used in this study.
Descriptive analysis. We used numerical and categorical "presence" of HD cases in a county. We excluded information from captive deer. We converted the numerical count data to "presence" and merged it with the categorical data from SCWDS , WHISPers (2004), and IDNR (2005) to produce an HD "presence" dataset. We screened for duplicate records across sources and removed duplicates.
Although HD data collection in Illinois began in 1982, the first reported cases of HD did not occur until 1988; therefore, we included only the years 1988 to 2019 in our analysis of HD presence.
We organized the Illinois counties by latitude from north to south using the county centroid and provided a visual representation of HD presence (Fig. 1) and numerical data (Fig. 2) for HD cases in wild WTD.
Disease mapping. Disease mapping, kriging, and the local hot spot analysis were performed using Spatial Statistics and Geostatistical Analyst Tools in ArcGIS 10.7.1 (Environmental Systems Research Institute, Inc., Redlands, CA, USA).
All reported HD cases during the study period were geocoded and spatially joined to the county level. The total number of HD cases by county was illustrated in ArcGIS 10.7.1 in a choropleth map.
Spatial interpolation of the county-level total number of reported HD cases was done using the Empirical Bayesian kriging method. To measure the spatial dependence of HD cases across Illinois counties, kriging uses a semivariogram, a function representing the distance and direction of two locations (e.g., counties). The Empirical Bayesian kriging accounts for the error created by estimating the semivariogram parameters by building several semivariograms instead of one and uses a restricted maximum likelihood estimation 76,77 .
Hot spot analysis of county-level HD cases was conducted using the Getis-Ord Gi* method 78 . Each county was represented by a polygon, its centroid, and its total number of reported HD cases for this analysis. The null hypothesis assumes complete spatial randomness of counties with high or low reported HD cases. The statistic provides a Z-score and p-value, which are associated with the standard normal distribution. High (positive) or low (negative) Z-score, associated with small p-values, are observed in the tails of the normal distribution and signify a statistically significant pattern. A significant (p ≤ 0.05) hotspot area is indicated by a high positive Z-score (≥ + 1.96) where counties with a high number of HD cases are surrounded by counties with a high number of HD cases. Likewise, a statistically significant cold spot (p ≤ 0.05) area is signified by a high negative Z-score (≤−1.96) and is detected if counties with a low number of HD cases are surrounded by counties with a low number of HD cases. To adjust for zoning (e.g., two counties might have the same number of cases but the county shapes and boundaries differ) and edge effects (e.g., counties at the periphery of the study area have neighboring areas that are not included in the study region), the "zone of indifference" conceptualization parameter was used.
The spatial interpolation and hot spot analysis results of county-level HD cases were illustrated in ArcGIS 10.7.1. The space-time permutation model requires only case data without a need for additional populationat-risk data. The model identifies disease clusters by comparing the observed number of cases in a specific area during a period to the number of expected cases that would occur if the spatial and temporal locations of all cases were independent of each other. It does not make prior assumptions on the size, location, or period of the outbreaks. It accounts for purely spatial and temporal variations in disease rates. For example, clusters will not be identified if during a specific period all regions have an increased number of cases than normal. However, clusters will be identified if cases increased in a specific region during a specific period and all the other areas have normal case numbers.
Space-time scan statistic. Retrospective space-time analysis scanning for clusters with higher than expected HD cases using the space-time permutation model 36 was conducted in the SaTScan software version 9.6 79 . The outcome variable was represented by the number of HD cases reported in a county in a year. The smallest spatial unit was signified by the geographical center (i.e., centroid) of a county, and the smallest time unit was signified by the year of the report. Cartesian latitude and longitude coordinates for each county's centroid were calculated in ArcGIS 10.7.1.

Scientific Reports
| (2022) 12:6888 | https://doi.org/10.1038/s41598-022-10694-y www.nature.com/scientificreports/ The space-time scan statistic uses a cylindrical scan window with a circular spatial base and height corresponding to time 80 . A scanning window containing up to 50% of the study area and up to 50% of the study period was used. The scanning window of variable radii gradually moves through time and space, evaluating the number of HD cases inside the scanning window comparing it to the number of cases outside the window. A primary cluster with the highest significant (p ≤ 0.05) test statistic after a Monte Carlo hypothesis testing with 999 replications is identified. Secondary clusters are reported if they are significant and do not overlap with the primary cluster.
All significant primary and secondary space-time HD high-rate clusters were mapped in ArcGIS 10.7.1 (ESRI Inc., Redlands, CA, US) at the county level.

Data availability
Data supporting this study are available in the Supplementary Dataset.